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Abstract 

For systems in an externally controllable time-dependent potential, the optimal protocol mini- 
mizes the mean work spent in a finite-time transition between two given equilibrium states. For 
overdamped dynamics which ignores inertia effects, the optimal protocol has been found to involve 
jumps of the control parameter at the beginning and end of the process. Including the inertia 
term, we show that this feature not only persists but that even delta peak-like changes of the 
control parameter at both boundaries make the process optimal. These results are obtained by 
analyzing two simple paradigmatic cases: First, a Brownian particle dragged by a harmonic optical 
trap through a viscous fluid and, second, a Brownian particle subject to an optical trap with time- 
dependent stiffness. These insights could be used to improve free energy calculations via either 
thermodynamic integration or "fast growth" methods using Jarzynski's equality. 

PACS numbers: 05.40.-a, 05.70.-a, 82.70.Dd 
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I. INTRODUCTION 



The free energy difference AF between two equilibrium states is an important quantity 
in isothermal statistical mechanics. Strategies to extract AF from experiments or computer 
simulations are traditionally based on either thermodynamic integration or thermodynamic 
perturbation 1 which use one infinitesimally slow transition or many infinitesimally fast tran- 
sitions, respectively, between the two equilibrium states. A decade ago, Jarzynski proposed 
the remarkable relation 



which interpolates between these extreme cases using nonequilibrium work values W ob- 
tained from trajectories of finite time transitions between the equilibrium states at temper- 
ature T (with Boltzmann's constant ks = 1 throughout the paper)^. This exact relation 
which holds for any time-dependent driving described by an external control parameter X(t) 
has been extended to various fluctuation theorems^ 5 -^^. Although these (necessarily irre- 
versible) finite time transitions occur in nonequilibrium, the equilibrium quantity AF can be 
inferred from a sufficient number of trajectories either from computer simulations^ 1 ^^ 2 - 1 ^ or 
real experiments^^ 5 -. However, the convergence of the involved exponential average causes 
problems for far out of equilibirium transitions where the work W is substantially larger 
than the free energy difference AF— . In this regime, the exponential average is dominated 
by low work values which are very rarely sampled 17 . As a remedy, several path sampling 
techniques biasing the dynamics for low work have been proposecr^^. It is, however, still 
under debate 2 ^ 2 - 1 ^ for which systems fast growth techniques are superior to refined "con- 
ventional" approaches such as umbrella sampling^ or flat histogram methods 25 . Though 
valuable for computer simulations, it is hard to imagine how to bias dynamics in real ex- 
periments, where, however, apparatus drift may prevent long measurements necessary for 
thermodynamic integration^ 5 -^ and thus render fast growth methods competitive. 

Both for thermodynamic integration and "fast growth" methods employing Jarzynski's 
equality (or some variant), efficiency gains can be achieved by optimizing the driving scheme 
A(i). For thermodynamic integration, where the work W > AF is taken as an estimator 
for AF, it is obvious that a minimal work gives the best result 27 . In the case of fast growth 
methods, the statistics for free energy estimates quite generally also improves with smaller 
mean work— Fast growth simulations even allow to combine data from different driving 
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schemes in a straightforward w&y—. 

The minimization of the work spent in a finite time process can, however, also be seen 
in the context of minimal energy dissipation. On a macroscopic scale, the optimization of 
the work (or power) exerted in a macroscopic cyclic process has been discussed for quite a 
while^* 3 ^ 3 - 1 - 13 ^. On a microscopic scale, fluctuations will also affect optimal cyclic processes 3 - 3 - 
which may become relevant for constructing optimal nano machines. 

For overdamped Langevin dynamics, the optimal protocol leading to a minimal mean 
work in a finite time t has been calculated analytically for harmonic potentials 34 . Surpris- 
ingly, the optimal protocol shows jumps at the beginning and end of the finite time transition. 
Since most molecular dynamics (MD) simulations of the dynamics of biomolecules are on 
time-scales where inertia plays an important role (see 3 ^ for a review on steered MD), it is 
an interesting question how these results transfer to underdamped dynamics. In particular, 
it is important to know whether the jumps are a result of having neglected inertia. 

In this paper, we investigate the role of inertia for two previously introduced paradigmatic 
processes. In Sect. [TH we calculate optimal driving schemes for an underdamped Brownian 
particle dragged through a viscous fluid by harmonic optical tweezers. In Sect. IIIIl we study 
an underdamped Brownian particle subject to an optical trap with time-dependent stiffness. 
In both cases, we compare our findings with the corresponding results in the overdamped 
limit 34 . We find that the optimal protocol still involves jumps. Even more surprisingly, the 
optimal protocol includes delta peaks at the beginning and end of the process. 

II. CASE STUDY I: THE MOVING TRAP 
A. Optimal protocol 

We consider a Brownian particle of mass m dragged through a viscous fluid with friction 
coefficient 7 by a harmonic potential 



where k is the (constant) trap stiffness. The focus of the optical trap A(t) is changed time- 
dependently from an initial position Aj = to a final position A/. Including inertial effects, 
the Langevin equation reads 




(2) 



mx = — 7X — k(x — A(t)) + rj{t), 



(3) 
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where thermal fluctuations are modeled by Gaussian white noise 



( V (t) V (t'))=2T 1 5(t-t / ). (4) 



The mean work spent in the process of total duration tf is given by 



^/ /dt (^M), (5) 



o 



where the average (. . .) is over the intitial thermal distribution and over the noise history. 
In the present case the total (mean) work reduces to 

W = ( f dtkX(X-u) = k { Xu+-X 2 f -k[Xut f , (6) 
Jo Jo 2 3 

where, for simplicity in the notation, we have defined the mean position of the particle as 
u = (x). This quantity u(t) depends on the whole history of X(t) and thus, the work W is 
a non-local functional of the protocol X(t). However, in analogy to the overdamped limit^, 
we can express the work as a local functional of the mean particle position u. By averaging 
the evolution equation ([3]) we have 

A = u + ju/k + mu/k, (7) 

which inserted in Eq. ([6]) leads to 



W 



m 2 ^ 2 2 m l I 2 2~l* 
2 U + 2k U + — UU+ 2k U 



J 



+ 7 / dtii 2 . (8) 



o 



The only term remaining in the integral, ii 2 , is identical to the one in the overdamped 
limit, while the boundary terms are different. In complete analogy to the overdamped case, 
we now proceed in two steps. First, we calculate the optimal shape u(t) minimizing only 
the integral given initial values u(0 + ) = and w(0 + ) = A. Note that despite the initial 
equilibrium value w(0) = 0, we are free to choose u(0 + ) = A since the necessary "kink" in 
u(t) at t — does not contribute to the integral. Similarly, at the end of the protocol (at 
t = tf) there can be another jump in the velocity. In a second step, we adjust the constant 
A to yield the minimal total work. First, from the Euler-Lagrange equation corresponding 
to the Lagrangian u 2 (and subject to the initial conditions just mentioned), we find 

u{t) = At (9) 
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for < t < tf. In contrast to the overdamped case, we cannot determine all the boundary 
terms at tf from the evolution equation. Thus, C = iiitf) is another free parameter. With 
u(tf) = [k(Xf — Atf) — 7(7] /m, we get the total work as a function of the yet unknown 
constants A and C 

W(A, C) = + \{\f~ Atff + 7 fj dtA\ (10) 

The work is clearly minimal for C* = 0, where the asterisk will denote optimal from now 
on. The remaining terms then read 

W(A) = t(\ f -At f y + 1 t f A\ (11) 

which, surprisingly, is exactly the same expression that was found in the overdamped limit. 
Minimizing this expression with respect to A leads to 

A* = — ^ (12) 

2 1 /k + t f y } 



which yields the work 



W* = kXl — -. (13) 

/ 2 + %/ 7 v ; 



Inserting Eq. (jUJ) into Eq. ((Tj), we find the optimal protocol 



for < t < tf implying symmetrical jumps 

AA = A(0+) — \ = Xf — X(tj) = Xf^j-^ (15) 

at the beginning and at the end of the process. 

Superficially, this optimal protocol looks like the expression in the overdamped case^. 
There is, however, a subtle difference arising from the presence of inertia terms. The optimal 
protocol forces the mean velocity to instantly jump at the beginning of the process from its 
initial equilibrium value u(0) = to ii(0 + ) = A*. At the end of the protocol, the optimal 
strategy consists in setting back the mean velocity to zero u(tf) = C* = 0. Due to the second 
time derivative in the equation of motion such jumps in the velocity, which require delta 
functions in the acceleration, imply a delta-type singularity in the protocol. Specifically, in 
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Eq. ([7]), the jumps in u imply a 5-function for u and hence a 5 function in X(t). The optimal 
protocol [Eq. ffT4|) ] thus becomes 

y(t) = \ f ptp± + -^f- [6(t) - 6{t - V)] , (16) 

as shown in Fig. HJ In the overdamped limit, m — > 0, the delta peaks vanish. 
B. Physical origin of singularities in the optimal protocol 

The benefit of having jumps in the optimal protocol can be understood intuitively as 
follows. From the perspective of minimal dissipation, it is obvious that the particle should 
be dragged at a constant (mean) velocity from the beginning rather than being accelerated 
during a finite time. This initial jump in the velocity of the particle can only be achieved 
by a finite initial difference A(0) — u(0), corresponding to a jump in A at t = 0. In the 
present underdamped regime, a velocity jump corresponds to a jump in the (mean) particle 
momentum which can only be achieved by a delta peak in the force, corresponding to a 
delta peak in the protocol. 

The final jump is harder to grasp intuitively. In fact, it stems from focussing on the 
minimal work rather than on the minimal (mean) dissipation (or entropy production). If we 
had searched for the minimal entropy production (as defined in 8 ), we would have found an 
optimal protocol without a final jump. In the present minimization, at the final time tf, the 
particle is still in non-equilibrium with respect to the final potential V(x, X(tf)). Relaxation 
to equilibrium leads to further dissipation after time tf which has, however, already been 
paid for by the total work since at constant A no work is exerted anymore. A smaller 
final particle position u(tf) leads to a longer relaxation time which can decrease the total 
dissipation of the combined process (nonequilibrium transition and relaxation). 

The final delta peak corresponds to setting the final velocity to zero. This decreases the 
kinetic energy of the particle and thus is beneficial for a small work. It also explains the 
surprising fact that, according to Eq. (113I) . we do not have to pay any extra cost for having 
inertia. During the initial singularity, the exerted work is stored in the (mean) kinetic energy 
of the particle. This contribution is fully recovered during the final singularity where the 
kinetic energy of the particle is set back to the equilibrium value. 
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C. Comparison to a linear protocol 



Without prior knowledge, one might have expected a continuous linear protocol 



A lin (t) = Xft/tf 



(17) 



to yield the lowest work. In the overdamped limit, the work for a linear protocol was at 
most 14% larger than for the optimal protocol. We now check how much smaller the value of 
the optimal work W* is compared to a linear protocol if we include inertia. First, we rescale 
the system in order to compactly write the relevant combination of parameters. With the 
rescaled mass m = mk/'-y 2 , the energy scale e = kX 2 and a rescaled time t = tfk/j, the 
work can be written as W = eW(t, rh), with the optimal work W* = e/ (2 + tj. 

Solving the second order differential equation of motion ([7]) using the linear protocol 
A lm (t), we find the ratio: 



W 



lin 



w* 



2+t 
P 

2+t 
i 2 



e 2m 

t_ 

e 2m 



9q cosh [utj + 6\ sinh(z4) 
6*o cos (ut) + 61 sin(z/t) 



m < 



m > 



(18) 



with 



4m — 11 



2m 



(19) 



and 



m 



1, e 1 



3m — 1 
2mu 



(20) 



In Fig. [21 we plot the ratio V^ lm /jy* as a function of rescaled time t and mass rh. This 
result shows that the optimal protocol significantly reduces the work spent in the process 
compared to a linear protocol. 



III. CASE STUDY II: THE STIFFENING TRAP 

In the first case study, only the averaged quantity u = (x) appeared in the work and thus 
the same result could have been obtained from a deterministic damped dynamics. We next 
examine a second case study where fluctuations are important. We consider a Brownian 
particle of mass m in an optical trap with time dependent stiffness X(t) which is driven from 
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an initial value A(0) = Aj to a final value X(tf) = A/ in a finite time tf. The time dependent 
potential 

V(x,t) = ^x 2 , (21) 
leads to the underdamped Langevin equations 

x = p/m (22) 
V — — r yp/m — \(t)x + r}(t), 

where p is the momentum of the particle and the noise i] (t) has the same properties intro- 
duced in the first case study. Again our main goal is to find the protocol for which the 
corresponding total (mean) work 

W = f tf dtX^- (23) 



o 2 

is minimal. Note that the mean squared position 

w = (x 2 ) (24) 
of the particle is non-trivially coupled to the mean squared momentum 

z = (p 2 ) (25) 
and to the position-momentum correlation 

y=(xp). (26) 

Their time evolution is governed by the set of coupled differential equations 

w = 2y/m, (27) 
z = -2Xy - 2-fz/m + 2 7 T, (28) 
y = z/m — Xw — -yy/m. (29) 

Unlike both the moving trap (with and without inertia) and the stiffening trap in the 
overdamped limit, the present case is much more involved since one cannot eliminate the 
protocol and write the work as a function of one variable only. We thus express the work as 
a time- local functional of w(t) and z(t). Solving Eqs. (1271) and ([291 for A yields 

A = —[z/m -jw/2 — mw/2] (30) 



which, inserted in Eq. (123]) and after partial integration, leads to 



W 

with the "Lagrangian" 



Aw 


mw 2 


tf 1 f 


- + 








8w _ 


o + 2^o 



dtC (31) 



•yw 2 zw mw 3 , . 

C = - + -r^r- 32 

2w mw 4w 2 



We proceed in two steps analogously to the moving trap. We first minimize the integral 
in Eq. (131]) for given initial conditions and then optimize with respect to remaining free 
parameters. The integrand £ depends on w (and w) but also on z. The variables w and 
z are not independent. Eliminating A and y from the equations of motion (!27j) . (128]) . and 
( |2~9]) . we find the physical constraint 

Q = zw — myw 2 /2 — m 2 ww /2 — 2yTw + 2'jwz/m + wz = 0. (33) 

A detailed analysis of the solution of this optimization problem using Euler-Lagrange equa- 
tions is given in Appendix IA1 

The results for both the rescaled optimal protocol \*(t/tf)/ Aj and the optimal work 
W*/T depend on the dimensionless quantities 

t = t f \i/j, X = Xf/Xi, m = Xim/-f 2 . (34) 

An extensive analysis of the optimal protocol as a function of all three parameters is out 
of scope. Since the overdamped limit (m —* 0) has been discussed previously, we focus 
on the behaviour as a function of m for given A = 2, t = 1. Given these parameters, the 
optimization problem can be solved numerically and the corresponding total work W can 
be calculated. In Fig. [3^, we plot the value of the minimal work W* (obtained from the 
optimal protocol) as a function of the rescaled dimensionless mass m and compare it to other 
benchmark protocols. All work values are bounded from below by the free energy difference 
AF = (In 2/2) T ~ 0.35T. Quite generally, work values are also bounded from above by the 
work for an immediate jump W^ p = lim 4 ^ W* = T/2. We study (i) a linear protocol, (ii) a 
protocol leading to a parabolic mean-squared position 

w (t) = ^(l + c t) 2 (35) 
Aj 

with optimized parameter c and optimized final delta peak, and (iii) a protocol leading to 

W(t) = l+at 3 (l- e -V[0.01+(5t/t,)»n + U 5 + cf 7 + df 9 (36) 



without any discontinuities (except for a final jump) but with free parameters a, b, c, d. The 
work arising from protocol (i) lies significantly above the optimal protocol. Protocol (ii) 
implies (optimized) jumps and delta peaks at the beginning and end. The work for protocol 
(ii) and the optimal work almost coincide. The inset shows that the optimal work is in 
fact slightly smaller than the work obtained for the protocol (ii). The difference to the 
numerically obtained exact solution W* decreases for decreasing m which is consistent with 
the analytical finding that protocol (ii) is optimal in the overdamped limit. The fact that 
protocol (ii) which involves optimized singularities but not the optimized shape is so close to 
the optimal work highlights the importance of jumps for the optimal protocol. Protocol (iii) 
has no delta peaks and no initial jump but mimics these features approximately since the 
parameters a, b, c, d have been optimized, see Fig. [3b. These trial protocols show that jumps 
and delta peak-like singularities can decrease the total work and confirms that our numerical 
solution of the Euler-Lagrange equations is the solution of the optimization problem. 

Finally, the explicit shape of the optimal protocol A* (t) can be reconstructed numerically 
from Eq. (1301) . see Fig. [3b. It displays initially a delta peak upwards accompanied with 
a jump AA and, finally, a delta peak downwards together with another jump AA'. Such 
discontinuities in the protocol are a consequence of the discontinuities in z, w and iv. The 
first singularity is "needed" to suddenly increase (p 2 ) from its equilibrium value and also to 
change the derivative of (x 2 ), which is proportional to the correlation (xp). Note that both 
size and direction of the jumps strongly depend on the rescaled mass fh as shown in Fig. 
[3H. For small m, the protocol jumps upwards (as also observed in the overdamped regime^) 
while for large rh, the protocol jumps downwards. 

IV. CONCLUDING DISCUSSION 

In summary, we have calculated optimal protocols yielding the minimal mean work for 
underdamped Langevin dynamics in two different model potentials. Surprisingly, these 
optimal protocols involve jumps and delta peaks at the initial and final times t j = and 
tf. While we have shown that the singularities in the optimal protocol appear for harmonic 
potential, there is no reason to believe that this feature generically vanishes for anharmonic 
potentials. In fact, in the overdamped limit, a recent study has shown that initial and final 
jumps are also present in a simple anharmonic potential 36 . At first sight, such singularities 
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seem to be unphysical since neither jumps nor delta peaks can be implemented in real 
experiments. Still our theoretical result is an important insight because it implies that 
there exists no optimal continuous protocol. Every such protocol could be improved by even 
steeper gradients mimicking the jumps and delta peaks at the beginning and end. If there 
was an experimental constraint on the allowed maximum rate of change in A, |A| < r, the 
minimal work would still be achieved by a protocol which looks roughly like the optimal 
one, with the jumps and delta peaks replaced by their best approximation consistent with 
| A | < r ( e. g. steep straight lines instead of jumps). Thus, it should be possible to exploit 
our results for real experiments. 

Our results may also be used in steered MD simulations. Even though we here have 
calculated optimal protocols for underdamped Langevin dynamics, there is no reason to 
believe that other thermostats frequently used in MD simulations would yield qualitatively 
different results for the optimal protocol. We have neglected memory effects by assuming 
white noise. While there are systems for which the underdamped Langevin equation is an 
appropriate physical description^, it might still be interesting to see how our results are 
altered when considering memory effects. 

The optimal protocol for a minimal mean work is not strictly equivalent to a protocol 
leading to an optimized free energy estimate. However, it has been found that the latter 
shares the same features (jumps at the boundaries) for overdamped Langevin dynamics^. 
Moreover, the optimal protocol leads to improved estimates of the free energy difference 
in both of our (underdamped) case studies. For case study I, the work distribution is 
Gaussian^ 1 ^ for which it has been shown that the error in the estimate of the free energy 
difference decreases with decreasing mean work—. In our second case study, for which the 
work distribution is no longer Gaussian, the error in the estimate of the free energy difference 
is indeed lower for the optimal protocol compared to a linear protocol, see Fig. H] and its 
caption for technical details. For thermodynamic integration, it is obvious that a minimum 
mean work leads to the best estimate of the free energy difference. We thus conjecture that 
appropriate singularities at the boundaries generically improve free energy calculations from 
either fast growth methods or thermodynamic integration. 

For determining the optimal protocol for an unknown potential, we envisage an adap- 
tive procedure in which trial protocols (including estimated singularities) are successively 
improved in an iterative fashion guided by the monitored work values. It might also prove 
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beneficial to use the optimal moving trap protocol (case study I) rather than a linear protocol 
in simulations of (protein) pulling experiments. 



Appendix A: SOLUTION OF THE OPTIMIZATION PROBLEM IN CASE 
STUDY II 

In this appendix, we give a detailed analysis of the numerical solution of the optimization 
problem. In order to minimize the integral in Eq. f[3"Tj) . the constraint [Eq. (|33)) ] is included 
in an effective Lagrangian £ e g- = £ — a(t)Q through a Lagrange multiplier a(t). Then the 
Euler-Lagrange equations whose solutions minimize the integral in Eq. (!3~T1) are obtained 
from 

dC c g ^ d 2 d£ c ff _ d dC c g dC e s _ d d£ c g (Al) 
dw dt 2 dw dt dw ' dz dt dz ' 

which, together with the constraint Q = 0, define a system of three differential equations for 
w, z and a. By defining the useful new variable 

2 

/i = zw —w 2 (A2) 

we can write the initially cumbersome differential equations ( 1A1I) after a tedious manipula- 
tion in the following reduced form 



w = — -- + 2Twa + — , (A3) 



w 2 2 a 2T 

-^ + 2Twa + — 

2w m z w m 

27 

ii = -u + 2-fTw, (A4) 

m 

2^i \ w ... 

a = — a H -. ( A5) 

m mw z 

These equations have no analytical solution but they can easily be solved numerically for 
given initial conditions w(0 + ), w(0 + ), /i(0 + ) and a(0 + ). It is important to note that some 
of these initial conditions are not fixed by the initial equilibrium conditions w(0) = T/Aj, 
w(0) = 0, w(0) = 0, z(0) = mT, but can be realized by additional discontinuities in the 
respective quantities at the boundaries. If such discontinuities do not change the value of 
the integral in Eq. fl3T|) . they do not affect the optimization of the integral via the Euler- 
Lagrange equations and hence the respective initial conditions should (in a first step) be 
treated as free parameters. Since the Langrangian does not depend on w(t), discontinuities 
in w{t) and w(t) can occur at the boundaries. However, a jump in the mean squared position 
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w(t) would affect the integral in Eq. (13"I1) and thus w(t) must be chosen to be continuous 
at the boundaries, enforcing w(0 + ) = w(0) = wq = T/Aj. Likewise, discontinuities in z{t) 
can occur at the boundaries. However, the initial values z(0 + ) and w(0 + ) are related by the 
constraint Q = 0. Integrating this constraint 



lim f t+e dt'G = (A6) 



leads to 



t + m" r nlt+ 

■ w 

J t- 



(A7) 



When applied at t — it yields 

T TO 2 

-[^(0+) - mT] = — w 2 {0 + ). (A8) 
Aj 4 

We consider a (possible) discontinuity through the parameter s\ in 

z(0 + ) = mTs x . (A9) 
With Eq. (IA8I) . the jump in the derivative of to at the initial time as a function of Si becomes 



w(0 + ) = ±27y^— i (A10) 
V mAj 

In the case in which Aj < A/, the correct sign is the negative one. Note that the last equation 
implies s\ > 1, so that at the initial time and given the equilibrium initial distribution, it is 
not possible to have a decrease in the mean squared momentum. From Eqs. (IA8I) and flA2j) 
we also find 

/i(0 + ) = mT 2 /\i. (All) 
Secondly, we define a new free parameter S2 in 

y(0+) = Ts 2 , (A12) 

which, from the evolution equation ( !29l) . directly yields iu(0 + ) = —S2- Then, writing Eq. 
( 1A3I) at t = + and inserting the above values, the initial value of the Lagrange multiplier 
needed to solve the Euler-Lagrange equations is 

a(0 + ) = A( aa _ ai + i). (A13) 
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Last, from the evolution equations (|27|) - (|29|) we find the relative value of the initial jump in 
the protocol as a function of si and S2: 

A(0+ 



A, 



s 2 - 1 + 71 



Isi - 1 



(A14) 



At the end of the process, the value of z(tf) is allowed to jump again. Recalling Eq. flA7[) 
applied now at the final time t = tf and isolating z(tf), we obtain 

m 2 w{t f ) 2 - w{r f f 



z(t f ) = z(t- f ) + 



(A15) 



4 wfo) 

Every quantity on the right hand side of the last equation except for w(tf) is fixed by the 
solution of the Euler-Lagrange equation. The minimum value for z(tf), which leads to the 
minimal contribution to the work in Eq. flA18j) . is reached for u>(i/) = s 3 = 0. 

For a comparison of the present case with its overdamped analogue^, one can formally 
integrate the differential equations for fi and a and plug them into Eq. (IA3I) to obtain the 
following integro-differential equation for w, 



w 



w 
2w 



IT 

m 



f(t)A--^ r) +f(t)(l + s 2 -s 1 ) 



where f(t) = ^ e 2 ^ m and 



A 



27 



-dt' B = 1 + 



27 



f(t')dt'. 



(A16) 



(A17) 



m Jo f(t') m Jo 

In the overdamped limit, the Euler-Lagrange equation is given by w — w 2 /2w = 0. Including 
inertia leads to nonvanishing terms on the right hand side of Eq. (1A16I) . However, taking the 
corresponding overdamped limit m — > in Eq. flA16[) yields the overdamped Euler-Lagrange 
equation only after optimizing the parameters s± and &2- 

Combining Eqs. (1271) -(|29l). the work W [Eq. ( l23i) ] can be written as 



W 



Xw z 



IT 



7 



tf H I dtz. 



711 



m 2 Jo 



(A18) 



To calculate the integral, we insert the solution of the Euler-Lagrange equations for z, which 
depends on si and S2- Then, we need to insert the boundary values for w and z at t = 
and t = tf. In a last step, the work is optimized with respect to the free parameters s\ and 
s 2 - 
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Figure captions 

Figured]: Scheme of the optimal mean position u*(t) and protocol X*(t). 



Figure [2] : Ratio between mean work W hn spent using the continuous linear protocol 
A lm (t) and optimal work W* as a function of the dimensionless parameters rh = mk/'j 2 and 

i=t f k/r 



Figure [3] : Optimization results for case study II for t — 1 and A = 2. (a) Mean work W* 
in units of T as a function of the rescaled mass rh compared to (i) a linear protocol, (ii) a 
protocol leading to w(t) given by Eq. (I3"5"j) . and (iii) a continuous (except for a final jump) 
protocol leading to w(t) given by Eq. fl36|) with adjusted parameters to yield a minimal 
work (see main text for details), (b) Protocol (iii) with optimized parameters for rh = 1. 
(c) Optimal protocol A*(t) for rh = 0.5. (d) Jump heights AA and amplitudes T> of delta 
peaks (in rescaled time tjtf) for the optimal protocol as a function of the rescaled mass rh. 



Figure H] : Comparison of free energy estimates for case study II for a linear protocol 
and a continous approximation to the optimal protocol for t — 1, rh — 1, A = 2. The 
data were obtained from Langevin simulation of 10 6 trajectories for each protocol with 
7 = 1, m = 1, T = 1, Aj = 1. (a) Distribution P(W) of work values W for the two pro- 
tocols shown in the inset: (i) the linear protocol A(t) = 1 + t and (ii) the linear proto- 
col with additional continuously approximated delta singularities. The columns show the 
free energy difference AF ~ 0.3466 and the mean work values (w®) ~ 0.4700(±0.0006), 
(WW) ~ 0.4270(±0.0005) which both are consistent with a direct evaluation based 
on Eq. (1231) . (b) Histogram of 10 5 Jarzynski estimates for the free energy difference 
AF est = — (l//3)ln J2iLi exp(— j3Wi)/N obtained from = 10 single trajectory work val- 
ues Wi each. The mean squared error (MSE) of these estimates consists of two parts^: the 
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systematic error (bias) B = (AF cst ) — AF and the statistical error a = yVar(AF cst ). The 
columns show the free energy difference and the mean value of the estimates obtained from 
the two protocols. Since the bias (B® = 0.0066(±0.0004), B™ = 0.0052(±0.0003)) can be 
neglected for both protocols, the MSE is dominated by the statistical error (a^ = 0.119, 
a {n) _ 0.104) which is smaller for protocol (ii). 
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